Spin-Orbit Resonance and the Evolution of Compact Binary Systems 
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Starting with a post-Newtonian description of compact binary systems, we derive a set of equations 
that describes the evolution of the orbital angular momentum and both spin vectors during inspired. 
We find regions of phase space that exhibit resonance behavior, characterized by small librations of 
the spin vectors around a fixed orientation. Due to the loss of energy and orbital angular momentum 
through radiation reaction, systems can eventually be captured into these resonance orientations. 
By investigating the long-term evolution of compact binaries with a variety of initial conditions, 
we find that the distribution in parameter space can be strongly affected by resonance captures. 
This has the effect of significantly reducing the size of search space for gravitational wave sources, 
in turn improving the chances of detecting such sources through methods of template matching. 
Furthermore, by calculating the expected spin distribution at the end of the inspiral phase, we can 
predict what are the most likely initial conditions for the plunge phase, a result of great interest for 
numerical relativity calculations. 

PACS numbers: 04.25.Nx, 04.30.Db, 95.30.Sf 



I. INTRODUCTION 



The inspiral and subsequent coalescence of two black holes (BH) or neutron stars (NS) promises to be a strong 
source of gravitational waves for a number of new interferometric detectors currently being developed 1 . In order 
to successfully detect and then analyze such sources, we need to have an accurate description of the gravitational 
| waveforms they will produce. For longer signals lasting many orbital cycles (typically NSs in the LIGO frequency 
band and BHs in the LISA band), the detection methods rely heavily on a technique called matched filtering. This 
ON ' method is based on the premise that we can calculate theoretical templates of gravitational waves which in turn are 
cross-correlated with the observed data in an attempt to match a small amplitude signal buried under high background 
noise. Since the anticipated signal can be hundreds or even thousands of cycles long, it is critical that the form of the 
theoretical template is very accurate or the detection could be missed. Once detected, the observed waveform must 
then be compared to a larger collection of model templates in order to fit the binary parameters and determine their 
Q-f statistical confidences. 

Approximate templates that do not include spin effects have been shown to have a poor chance of matching 
gravitational waveforms from spinning binaries Q, 0, 0, 0, B El • Even if we were able to calculate the theoretical 
waveforms with perfect physical accuracy, the binary black hole system is so complicated that we would need a very 
large template library to give a reasonable chance at detection 0,0. For two spinning black holes, the parameter 
space is characterized by at least 11 intrinsic variables [the masses (2), the angular momentum vector (3), and two 
spin vectors (6)]. Apostolatos et al. Q showed that the two-spin system can be reduced in the limits of equal mass 
when neglecting spin-spin interactions and later Apostotalos |9( included these terms for the equal mass, equal spin 
case. Recent work by Buonanno et al. @] showed that the size of this parameter space can also be reduced by 
considering a set of quasi-physical templates that mimic the physical behavior of two spins with a single effective 
spin. Grandclement et al. Q have suggested a different method of using "spiked" templates to expand the search 
templates and simulate spin effects. While these fitting methods greatly aid the searches for gravitational waves, they 
could also benefit from additional astrophysical information about the systems that are producing the waves, as well 
as their evolution up to the point where they enter the frequency regime of the detector. The added advantage of 
using strictly physical templates is the direct manner in which they allow us to determine the intrinsic parameters of 
the compact binary. 

There are a number of stellar evolution models that describe the formation of binary black hole systems, inclu ding 
estimates for initial spins and kick velocities, which in turn can give the orientation of the orbit 0, 0, IT2I Il3| . 
However, there is still a fair amount of uncertainty in the appropriate initial values to use for inspiraling stellar mass 



1 http://www.ligo.caltech.edu/ 
http:/ /www. virgo.infn.it 
http: / /www. geo600.uni-hannover.de 
http: / /tamago. mtk.nao.ac.jp 
http: / /lisa. nasa.gov 
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black hole binaries. The mechanisms governing supermassive or intermediate mass black hole mergers are even less 
certain. In both cases, we have little or no idea what to expect the system might look like as it enters the final stages 
of evolution towards inspiral and merger. 

In addition to the compact objects formed through binary stellar evolution, another important source of gravitational 
waves may be capture binaries |14| . These include binary systems that form in the cores of dense globular clusters, 
stellar mass and supermassive black holes in the centers of galaxies, and supermassive black hole binaries formed 
through hierarchical galactic mergers. Many of these systems should have no a priori preference for any particular 
spin orientation. It is also possible that some of these systems have already been seen, but not through the detection 
of gravitational waves. A recently discovered binary pulsar system ma y gi ve important information about the initial 
spin-orbit orientations and early evolution of inspiraling neutron stars |l5l fill . Merritt & Ekers ^3] suggest that the 
morphology of radio jets from active galaxies may point to recent supermassive BH mergers and resulting changes in 
spin orientation. In fact, the longer term merger history of a supermassive black hole might also be inferred by its 
current mass and spin UM- If we could give a realistic prediction of the orientation of both spins and the angular 
momentum vector relative to each other at the time of inspiral, an important link between these different branches 
of astrophysics could be established. A successful identification with an electromagnetic source would also be critical 
in the confirmation of any gravitational wave detection, especially in the early years of gravitational wave astronomy. 

To further investigate these issues, we have developed a code to integrate the post-Newtonian equations of motion 
and spin precession for two spinning point masses including radiation reaction. The full evolution of the system 
from its formation until merger covers an enormous range of time scales, from hundreds of megayears to fractions of 
milliseconds. To accommodate this range, we first derive an orbit-averaged system of evolution equations that model 
the orientation of the spins and angular momentum without actually following the orbital phase of the binary. These 
orbit-averaged equations of motion agree quite well with the full 2.5-order post-Newtonian equations for modeling the 
evolution of the orbital angular momentum and both spin vectors |19( . 

With this orbit-averaged formulation, a collection of equilibrium solutions are found in which the relative orientation 
of the spin vectors and orbital angular momentum remains fixed in time. Furthermore, we find that the majority 
of these equilibrium solutions are stable, so systems nearby in parameter space librate around the stable orbits like 
a spin-orbit or spin-spin resonance. With the inclusion of radiation reaction, spin-locked systems remain locked, 
following a trajectory along equilibrium solutions with steadily decreasing orbital angular momentum. Furthermore, 
initially non-resonant systems can actually be captured into these stable orbits and then oscillate around a fixed 
orientation throughout the rest of the inspiral process. 

One effect of this resonance behavior is to significantly reduce the size of the parameter space over which the 
spin orientations are distributed. In turn, this could have an important impact on the size of the template library 
used for waveform matching. Instead of including a separate template for every single spin orientation (at least four 
parameters), we could instead use a one- or two-dimensional family of equilibrium solutions as representative of the 
entire sample. These "guiding center" waveforms would be a good match for any resonant system, librating around the 
equilibrium orientations, potentially reducing the size of the search space by orders of magnitude while maintaining 
a set of strictly physical templates. 

Alternatively, we can take the inverse approach so that, given the gravitational wave form of an inspiraling binary, 
we should be able to infer the evolutionary history leading up to the point where it enters the detector frequency 
band. In this manner, perturbative methods could be used to explore the phase space around the equilibrium regions 
and better pinpoint the exact binary parameters of the detected signal. Eventually, as the detection rate increases 
to hundreds or even thousands of signals per year, we will be able to fully map out the spin distribution of compact 
binaries in the Universe and thus test the theoretical results presented in this paper. 

In addition to reducing the volume of parameter space in the template library, the equilibrium solutions could be 
extremely important for determining initial conditions for the plunge calculations. This phase of the binary merger is 
widely regarded as the least-well understood physically and at the present time we have relatively little faith in any 
theoretical waveforms that might describe the plunge (see, e. g. [2(| and references therein). Since the post-Newtonian 
expansion (known at this point up to terms of order 3.5PN) becomes less and less accurate in the plunging region, it is 
likely that we will need to rely on full numerical relativity calculations to produce the transitional waveforms linking 
the long inspiral and the ring-down phase (which is itself well described by perturbative methods). Because of the 
prohibitive computational expense of these calculations, it is that much more important to have a solid understanding 
of how to describe the initial conditions. We believe that the equilibrium solutions described in this paper are an 
excellent place to begin, as they may represent a large segment of the relevant binary population. Knowing the final 
spin orientation at the plunge phase would also improve our understanding of black hole recoil velocities, a field of 
much recent interest |21| . 

In Section [H] we present the evolution equations and discuss the relevant time scales for the binary inspiral and 
spin precession. In Section IIIII we show how this system of equations can be further reduced to five variables (four 
when assuming circular orbits), which completely describe the relative orientation of the compact objects, and how 
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this reduced system exhibits resonance behavior. Section llVl describes the evolution of such a system under radia- 
tion reaction and the dependence on initial conditions. In Sections IV! and fVTl we discuss the effects of large mass 
ratios, eccentric orbits, and the dependence on spin magnitude. We conclude with a discussion of the astrophysical 
implications of these results and directions for future work. 



II. BINARY EVOLUTION EQUATIONS 



We use the post-Newtonian equations of motion for the reduced two-body problem including Lense-Thirring pre- 
cession terms for the point mass spins. Throughout t his p aper we adopt geometrized units with G = c = 1. The 
instantaneous precession equations are (see, e. g. |& ll9ll22l |) 

Si = «i x Si (2.1a) 



where 



s 2 = n 2 x S 2 , 



(2.1b) 



ill 



3 m 2 
2 mi 



L w -S 2 + 3(f -S 2 )f 



(2.2a) 



fl 2 = ^7 



^)L Ar -Si+3(f -Si)f 

2 m 2 



(2.2b) 



Here Si i2 are the spin vectors of the two compact objects with units of angular momentum and a dot represents 
time derivative (assuming a global coordinate time t in the post-Newtonian approximation) . The bodies have masses 
mi > m 2 and are separated by a distance r in the f direction. The Newtonian orbital angular momentum Ljv is 
defined in the usual way: Ljv = ^t(r X r) with the reduced mass defined as 



mim 2 



(2.3) 



m 2 . 



and the total mass is m = mi 

At first glance, the leading r~ 3 terms in equations H2.2al2.2b|) suggest that spin precession only becomes important 
late in the inspiral when r becomes small. However, since the great majority of the evolution (in terms of time) occurs 
at a large separation, the spins actually have a chance to undergo a large number of slow precessions leading up to 
the inspiral phase. To first order, the evolution of a circular orbit due to gravitational radiation is |23| 



where a = ^-[im 2 , giving 



r(t) = (r 4 (0)-ai) 1/4 , 



We can thus define a characteristic time scale for inspiral as 



T, 



inspiral 



The precession time scale is given to leading order as 



dr/dt 



(2.4) 



(2.5) 



(2.6) 



(2.7) 



For large r, the precession time scale is actually shorter than the inspiral time scale, and thus we see that the spin-orbit 
effects are important even at early times. 
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At these large separations, the system evolves very slowly in time, with an orbital period much shorter than the 
precession period. The number of orbits iV or b at a given separation r and orbital period T or b is given by 



1 



dr 



T Q1 



while the number of precession cycles N pTec is given by 



dr 



T 

- 1 - pre 



r 3/2 



,1/2 



(2.8) 



(2.9) 



In short, this means we must calculate the spin evolution over a large range of r, corresponding to a very large number 
of orbits. However, since we are interested here in the net evolution of the spin orientation, the actual phase of the 
binary orbit can be ignored and equations (|2.2a(l and 12.2b|) can be calculated in an orbit-averaged approximation 
(see Appendix^!. Generalizing to non-circular orbits with semimajor axis a and eccentricity e, the orbit-averaged 
precession vectors are 



and 



1 



a 3(l_ e 2)3/2 



a 3(l_ e 2)3/2 



3 m 2 
2 mi 
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2 + -^1 - 3 Sl ' Ljv 
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1 
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(2.10a) 



(2.10b) 



An interesting result is that the orbit-averaged precession equations are independent of the ellipse's orientation in 
the plane normal to L^r. In other words, the first-order post-Newtonian (1PN) effect of pericenter precession does 
not play a role in the spin evolution, which is fortunate because of the relatively short time scales involved, compared 
to the slower effects of spin-orbit acceleration (1.5PN), spin-spin acceleration (2PN), and radiation reaction (2.5PN) 
0,122. In the orbit-averaged equations of motion, the spin-orbit precession terms (1PN) dominate on a short time 
scale, but the secular effects of the radiation reaction are also very important in the long-term evolution of the system. 

On time scales short compared to the inspiral time, the total angular momentum 



J — Ltv + S = Ljv + Si + S2 
is conserved, implying that the orbital angular momentum evolves according to 

Lat = — S ~ — S c ff x Ljv, 



(2.11) 



(2.12) 



where S e g is some linear combination of Si and S2 (note this is not the same effective spin as in 0). Since Ljv = 0, 
without radiation reaction, the magnitude of the orbital angular momentum vector is constant with spin precession. 
The magnitude of the total spin vector S, on the other hand, is not conserved as the angle changes between the two 
spin vectors (each of constant magnitude). 

The relations (|2.11|l and l|2.12|l constrain the binary system to a subset of the complete parameter space defined 
by the three vectors L^r, Si, and S2. We believe it is this set of constraints that best explains much of the behavior 
presented below, as opposed to a more classical description of resonance based on Hamiltonian mechanics and energy 
minima in phase space (see, e. g. Murray & Dermott [24|], Sussman & Wisdom j2^|). However, a Hamiltonian 
formulation of the post-Newtonian equations of motion such as in [26l I27I |28| may prove to give a more classical 
explanation to these apparently geometric constraints. 

The inclusion of gravitational radiation causes the orbit to shrink and also circularize in time, reducing a, e, and 
the magnitude of the angular momentum 



Ln = \i\J ma(l — e 2 ). 

Following Peters |23| and adopting units with m = 1, we use the coupled first order differential equations 



G4 



dt 



5 a 3 (l -e 2 ) 7 /2 
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(2.13) 



(2.14) 
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"iy fl 4 (1 _ e2) 5/2 ( v 1 + 5o4 e )' ( 2 - 15 ) 



<f T _ 32 m 2 ( A , 7 2 



and 

to evolve the binary orbital elements in time. All of the above orbit-averaged precession and radiation reaction 
equations have been tested and compared to the full 2.5-order post-Newtonian equations of motion in Kidder |19| . 
The agreement is very good for most of the inspiral, all the way down to r ^ 10m, after which almost any post- 
Newtonian approximation becomes increasingly uncertain. 



III. GEOMETRY OF EQUILIBRIUM 



One of the most difficult aspects of studying the spinning binary system is the problem of visualizing and analyzing 
the orientation of the two spins and the angular momentum in an informative way. In general, these three vectors 
are defined by nine coordinates [the angular momentum is also related to a and e through Q2.13[) ]. Since the spin 
magnitudes Si and S2 are conserved in the point mass approximation, and we can pick a coordinate system where L^r 
points in the e z direction, we are left with five coordinates: {L^, 9±, 0%, (pi, </> 2 ). Furthermore, the overall dynamics 
are preserved under rotation around Ljv so we can reduce the spin degrees of freedom by defining the e x direction 
along <j>i — 0, leaving four independent coordinates to define the orientation of the system: (Zjv, 61, Q%, A</>). Figure 
shows a schematic of the geometry used throughout this paper. Following the post-Newtonian formalism, all angles 
and vector magnitudes are defined in a Cartesian, flat space-time. 

In this coordinate system, there exists a set of equilibrium spin configurations for which Ln,9i,02, and A<p are 
constant (without radiation reaction), even though the individual vectors might vary in time from the point of view 
of a fixed inertial coordinate system. Trivial equilibrium examples include the collinear cases with cos#i = ±1 and 
cos 82 — ±1. More interesting cases occur when Si, S2, and Ljv all appear to precess around a fixed axis at a constant 
rate so as to remain in a fixed relative orientation. These points in parameter space can be found by solving (cf. 
Apostotalos 0) 



di (Sl ' S2) 2a 3 (l-e 2 ) 3 / 2 



mi _ m± 4. ( s i - s 2) • Lw 

mi m 2 L 2 N 



S 2 • (Ljv x Si) = 0. 



Note that the last term (a vector triple-product) can be written in our reduced coordinate system as 

S2 • (Lat x Si) = S1S2LN sin 6>i sin 62 sin A</>. 



(3.1) 



(3.2) 



Thus l|3.1|l is satisfied when Si, S2, and Ljv are coplanar, i.e. sin A0 = for all times. In practice this means finding 
simultaneous solutions to 



S 2 • (L w x Si) =0 (3.3) 

and 

^[S 2 -(L Jv xSi)] = 0. (3.4) 

Combining with equations (|2.2al l2~2b1 12.10al and l2.10Tijl . 13.4fl can be written in terms of just Si, S 2 , and Ljv, with 
no explicit time derivatives: 

(fii x Si) • [S a x (L N + Si)] = (n 2 x S 2 ) • [Si x (L N + S 2 )], (3.5) 

which in turn are simple combinations of our reduced coordinates (Ljv, &i,02, A0). 

For a given value of Ln (i- e. at a particular point in time during the binary inspiral), and setting sinA0 = 0, 
solutions of (|3.5|l trace out one-dimensional curves in (6 1,# 2 ) space. Figure [3 shows these curves for two maximally 
spinning black holes with nearly equal masses (mi = 0.55, m 2 = 0.45) at a few stages of the inspiral evolution 
(Ljv = 4,3,2, 1,0.5, with binary separations of r/m w 260, 150,65, 16,4). These solutions correspond to the aligned 
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configuration with A<fi — and are found numerically to be stable, as will be explained below. There also exist 
anti- aligned equilibrium curves for A</> = 180°, shown in Figure for the same masses and values of Ljv- We find that 
the curves along the bottom (0 2 50°) and the right side (0\ > 140°) of the plot correspond to stable equilibriums, 
while the solutions in the upper left corner, appearing at later times, are generally unstable or quasi-stable. 

The regions of stability can be better understood by considering a type of effective potential for orbits in parameter 
space near the equilibrium solutions. While the precession equations do not suggest any obvious way of being 
reformulated in the traditional approach of an analytic effective potential, we can map out the potential numerically 
by integrating a large sample of initial conditions. For each of these initial conditions, the stability of the spin 
orientation can be measured by integrating the change in A<p: 



With this definition, [A0(t)] to t is potentially unbounded and can grow in magnitude larger than 360° as one spin 
precesses around Ljv more rapidly and repeatedly "passes" the other spin vector. For binaries locked in a spin 
equilibrium (with either A<fi — 0° or A<f> — 180°), the integrated phase shift will be zero. For systems far from this 
equilibrium, the two spins will precess more or less independently and will generally acquire a linear phase shift in 
time since ^ And for systems near the equilibrium points in phase space, we expect the phase to librate 

around the exact coplanar solution given by equation (|3.5I) . 

Figure 0] shows [A</>(i)] to t for a few examples of these different regions in parameter space, again with rti\ — 
0.55, m 2 — 0.45, and maximum spin parameters (Sj = rrij). Using Figures |21 and [3] as guides, we can identify initial 
conditions near stable and unstable equilibrium solutions. Figure shows a binary near a stable equilibrium with 
9i = 50°, 6> 2 = 120°, A(f> = 0, and L N = 4. The system librates around A<j> = with amplitude of about 11°. For 
initial conditions far away from equilibrium points, the two spins precess at independent rates, giving a roughly linear 
drift in phase, as seen in Figure^, which has initial conditions 9i = 120°, f? 2 = 50°, A</> = 0, and Lm = 4. These two 
examples are also labeled on Figure|21as 'A' and 'B'. 

For the systems oscillating around a stable equilibrium, we can think of an effective potential whose form can be 
inferred by the libration of trajectories through those points in phase space. With this approach, the minima of the 
potential lie along the solid curves plotted in Figures [21 and |3J Systems lying near these curves in parameter space 
will oscillate around them with constant amplitude and frequency. Phase space trajectories like the one in Figure 2b 
correspond to a plateau-type region of the effective potential with constant amplitude, where the spin vectors can 
precess freely through all values of A<p. 

From this analysis, we can also determine that some of the equilibrium solutions plotted in Figure |21 are unstable 
(dashed lines). Perhaps it is more accurate to call these regions quasi-stable since initial conditions a small distance 
outside of these closed curves are found to exhibit behavior typical of trajectories near an unstable equilibrium: long 
periods of stasis followed by short bursts of rapid divergence away from equilibrium, as in Figure (9i = 30° , 6> 2 = 
133°, A(j> = 180°, Ljv = 1, labeled 'C on Figure|3J). However, points just inside the equilibrium curves appear stable, 
librating around the locked A(f> = 180° position in phase space. 

In practice, we find that binaries initially outside of these quasi-stable regions tend to remain outside. Furthermore, 
when including radiation reaction, all binaries start outside of these regions, which only begin to appear for relatively 
small values of Ln late in the evolution. Thus it would be quite unlikely for an astrophysically realistic binary system 
to be found on the stable side of any of the quasi-stable equilibrium curves. Yet these regions are still of significant 
physical interest, as we believe the boundaries between stable and unstable regions in parameter space may be relevant 
to the ongoing question of chaos in spinning compact binaries and the existence of positive Lyapunov exponents in 



In Section IlIII we assumed in the analysis that the orbital angular momentum Ljv was constant in magnitude, i. e. 
there was no radiation reaction damping the system through gravitational wave emission. Of course, in all physically 
realistic black hole binaries, gravitational radiation will play a major role in the secular evolution of the orbit. From 
the post-Newtonian analysis in Section [HI we see that both radiation reaction and spin precession are dynamically 
important through most of the long inspiral process. 

With the inclusion of radiation reaction, the angular momentum evolves according to equations (|2.12|) and i|2.16[l . 
changing direction due to Lense-Thirring precession and losing magnitude due to gravitational losses. The individual 
spin vectors still maintain constant magnitude under radiation reaction, but have access to a broader range of preces- 
sion angles since the total angular momentum J is no longer conserved. In this sense, radiation reaction removes an 




(3.6) 



such systems (see, e. g. p^.!^ 1 ). 



IV. SPIN EVOLUTION 
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integral of motion from the precession equations, allowing a greater region of phase space to be sampled throughout 
the inspiral evolution. 

For the stable equilibrium solutions to equation Ij3.5|l . we find that the binary systems that lie on or sufficiently 
close to one of the curves in Figures [21 and [3] will evolve along successive curves of decreasing Lm as the system loses 
angular momentum to radiation reaction. At the same time, these stability regions grow stronger, confining a wider 
range of spins to the equilibrium curves. Even for some initial conditions that begin far away from the stable solutions, 
precessing freely through a large range of A(f>, the effect of radiation reaction causes the spins to become locked and 
then librate around equilibrium. These systems undergo transitions from the "drifting" behavior of Figure 2t> to the 
locked condition of0^. Once the system becomes locked, we find it will stay locked throughout the rest of the inspiral. 

The net effect of this spin evolution can be seen clearly in Figures and We show sinusoidal projections of the 
spin coordinates (# 12 ,A(/>) for an initially random, uniform distribution of spin vectors S2- The polar angle (9 12 is 
defined as the angle between the two spins: #12 = cos _1 (Si • S2) and the azimuthal angle A<fi is defined as above. 
A sinusoidal projection maps a point with spherical coordinates (0, </>) onto a plane with x-y coordinates given by 
(</>sin(9, 9) so that equal areas in the plane correspond to equal solid angles on the sphere. This allows an accurate 
way of estimating the density of the distribution function in parameter space. 

In Figure[5]we take a sample of initial conditions with r/m = 1000, 9\ — 10°, and masses mj = 0.55, m 2 = 0.45. 
The relative spin orientation (#12, A(j>) for the smaller black hole has a uniform random distribution, as shown in the 
upper-left frame of Figure [SJ As the binaries evolve under radiation reaction, more of the initially unlocked spins 
become locked and then librate around A<p = 0, as seen in the second frame of Figure El corresponding to r — 250to. 
Once locked into an equilibrium spin orientation, the systems move along the curves shown in Figure |2 approaching 
the line ^ = ^1 as L« -» 0. This evolution is clearly evident in the last frame of Figure [5] as #12, A</> — > and the 
two spin vectors become closely aligned with each other (although not necessarily aligned with the orbital angular 
momentum). This result may have a important impact on the question of initial conditions for compact binaries 
entering the plunge and subsequent merger phase. While the general problem of two spinning black holes merging is 
still an open question in numerical relativity, it would certainly be a significant simplification if we can assume that 
the spins are initially aligned. 

If we take the initial orientation of the more massive black hole spin to be retrograde with respect to the orbital 
angular momentum (#1 = 170°) and again evolve a uniform distribution of S2 orientations, we find a qualitatively 
different behavior. Instead of getting locked into the A</> = equilibrium solutions of Figure [5] and evolving towards 
9\2 = 0, the binaries approach the A<fr — 180° solutions of Figure where there is no strong correlation between 61 
and 62 at late times. Thus as Ljv — > and A<j) — > 180°, the distribution in 612 is roughly uniform, with a trend towards 
anti-alignment with 9\2 > 90°. This is seen in the evolution of the spin distribution shown in Figure |SJ Similar to 
Figure [3 a sample of spin orientations (#12, A</>) with uniform distribution is taken for the set of initial conditions 
with r/m = 1000 and 6\ = 170°, plotted in a sinusoidal projection in the upper-left panel. As the systems evolve 
under radiation reaction, they approach an anti-aligned spin configuration with A<fr — > ±180° and #12 ~ 90° — 180°. 

We see from Figures and that the final spin configuration is strongly dependent on the initial orientation of the 
larger spin with respect to the orbital angular momentum. Figures and [3] serve to quantify these effects for different 
initial values of 9\. When 6*1 (to) ~ 0°, the final spins tend to be parallel, regardless of the initial orientation of S2, as 
shown in Figure[7] However, if the smaller black hole initially has a retrograde orbit, the final spins tend towards the 
anti-aligned orientation described above (Figure|SJ). For the intermediate cases with #i(£o) ~ 90°, there seems to be 
little net evolution in the relative spin distributions. 

Given an initial probability distribution function we can predict the final distribution of (#12, A<fi) by weighting 

the ensemble of distributions in Figures [7| and |H1 by f{0\). This initial spin distribution function has been the focus of 
much recent work in compact binary systems |10L 1 1 01 ] . Much of this work is based on the evolution of two large main- 
sequence stars that both eventually form compact objects through supernovae |3l|]. The resulting spin and angular 
momentum vectors are largely determined by the random kick given to each star during the asymmetric supernova 
explosions. Many of these kicks will disrupt the system entirely, while the undisrupted systems may end up with 
very different orientations than they had before the supernova. Using these binary evolution codes to give the initial 
angular momentum and spin orientations, we can then use the above results to evolve the compact binary system 
under radiation reaction into the detector frequency regime. According to Kalogera 01 1 a supernova kick velocity of 
Vk ~ 200 km/s should favorably produce a system with 6\ < 30°. In this case, it follows that the A0 = 0, #12 = 
equilibrium solution would well describe the orientation as the system enters the detector's frequency range. However, 
even for systems locked in equilibrium, the total spin vector may not be parallel to the angular momentum vector, 
and thus simple precession Apostolatos et al. Q of the orbital plane will still take place during the inspiral and must 
still be accounted for in the template libraries. 
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V. LARGE MASS RATIOS 



Preliminary results with large mass ratios (mi : m-2 > 1000) show a similar behavior as the near-equal mass 
systems, although only for cases where 9i(t ) is close to zero [or alternatively 6i(t a ) ~ 180°]. This might correspond 
to a supermassive black hole in the center of a galaxy surrounded by a large "accretion disk" of stellar-mass black holes 
and neutron stars 32]. If the smaller black holes begin with randomly oriented spins, then as they evolve towards 
merger, we find that 0i2 remains roughly constant but Acf> — > 0, as in the previous section. However, it is more 
likely that any massive stars formed in a self-gravitating accretion disk would not have random spin orientations, but 
rather be aligned with the angular momentum of the surrounding gas, and thus the central supermassive black hole, 
simplifying the problem even further. 

If, on the other hand, the stellar-mass black holes originate from an isotropic population in the galactic bulge 
and sink towards the central black hole via dynamic friction, they should have initial velocities (i. e. orbital angular 
momentum) and spins with uncorrelated, uniform random distributions. In this case, the alignment effects described 
above are unlikely to be important, and the final spin and angular momentum orientations will likely have a similarly 
random distribution at the time of merger. While the large mass ratio limits the effect of spin locking, at the same 
time the spin ratio scales as mfj/m 2 , so to a large degree the smaller black hole can be treated as a point mass moving 
in a stationary spacetime, greatly reducing the parameter search space in any case 0. 

For the moderate mass ratios (mi : ra-i « 1 — 10) for which resonance behavior will be important, one can estimate 
roughly the range over which the spin locking occurs during the binary inspiral by examining equation (|3.1|1 . To 
leading order, we can write the equilibrium condition as 

d ( m 2 mi Si cos 0i - S 2 cos 2 N n 

— (cos i2 ) oc h «0. (5.1) 

at \mi m 2 J 

Taking maximal spins and Ljy = /Urn 1 ' 2 / -1 ' 2 , we can solve for the separation at which the resonance behavior begins 
to dominate the dynamics [maximizing r over the solution space to (|5.1|) ]: 

. /Si cos 0i - S 2 cos0 2 \ 2 (m\ + ml\ 2 

riock/m « - 2 ~ — 2 ■ ( 5 - 2 ) 

\ mj — m-2 / \ mf — J 

Thus we see that for large mass ratios, the spin locking does not occur until very late in the inspiral, while for equal 
mass systems it is important even for a relatively wide separation. This result is consistent with the recent findings 
by Buonanno et al. p] that find the best fitting factors in the cases of nearly equal masses, further emphasizing the 
similarity between their effective spin method and the equilibrium solutions of this paper. 

The dependence of spin locking on the mass ratio is shown in Figure El for maximally spinning black holes with 
0i(*o) = 10°. For mi : m.2 = 0.55 : 0.45, the spin locking is very strong at the end of the inspiral (r/m = 10), as 
we saw in Figure |SJ With increasing mass ratios of mi : rri2 — 3 : 2, 3 : 1, and 9 : 1, the resonance behavior grows 
increasingly weaker, showing a smaller influence on the initially uniform distribution of (0i2, However, it should 
be noted that since the inspiral time scales as the higher mass ratio binaries also spend a longer time at smaller 
separations, producing many more cycles in the waveform at late times. Thus the equilibrium solutions may still give 
good estimates for waveform templates even for the moderate mass ratios that would be expected for NS-BH binaries. 



VI. EFFECTS OF ECCENTRICITY AND SPIN MAGNITUDE 



For the sake of simplicity, most of the computations carried out so far have had circular orbits throughout the 
inspiral. However, it is relatively straightforward to include a non-zero eccentricity and evolve it along with the other 
system variables according to equations Q2.14I 12.151 12.16(1 . As shown in Peters [23] , a high eccentricity has the effect 
of speeding up the inspiral time scale. Yet as we see from equations (|2.10al 12 . 1 b|) . eccentricity also increases the 
precession rate, accelerating the rate of spin evolution. 

For the nearly equal mass cases discussed in Section IIVI we find that including a high initial eccentricity (ao = 
2000m, eo = 0.9) reduces the overall effect of the spin locking, but only slightly. A sinusoidal plot of (0i 2 ,A<^) 
analogous to Figure shows the same qualitative behavior for eccentric orbits, with marginally more scatter in the 
final distribution. This suggests that the increase in the inspiral rate slightly prevails over the increased precession 
rate and the spins do not have as much time to get locked into their equilibrium states. 

While we have not yet discovered a BH-BH binary, the known NS-NS systems show a range of eccentricities up to 
e w 0.6, presumably due to the variation of initial kick velocities [l^. By integrating l|2.14|l and (|2.15fl . it is clear that 
all these known systems will circularize long before the spin locking becomes important, and certainly before entering 
the detector band [ll"|. 
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There is also the possibility of forming high eccentricity BH-BH binaries through many-body exchange interactions 
in the cores of dense globular clusters, which are then ejected and merge within a Hubble time ^4[. Yet even for 
these systems, the merger time scale is long enough that the circular orbit approximation is most likely a reasonable 
one. In the case of the high mass ratios described in Section high eccentricities are much more likely. However, 
as was mentioned above, in such a system the contribution of the smaller spin is negligible and the solar mass black 
hole can be treated accurately as a point mass. Furthermore, from 1(5. 2|1 we see that spin effects are important at 
a much later stage for high mass ratios, giving more time for the system to circularize, again justifying the circular 
orbit approximation. 

We have also investigated the effects of non-maximal spin parameters on the evolution of the spin distributions. 
While equation (|5.1|) suggests that the resonant locking effects might fall off for smaller values of Si and S2 , in practice 
we find a relatively weak dependence on the magnitude of the spin parameter for a/m — S/m 2 > 0.5. Figure ITUI 
shows the probability distributions of (#12, A</>) near the end of inspiral for a variety of spin magnitudes. Each system 
has the same mass ratio mi : m,2 — 11 : 9 and initial spin-orbit angle #i(io) = 10°. Even for spins as small as 
a/m = 0.25 we see significant spin locking at the end of inspiral. Interestingly, the final distribution for a/m = 0.5 
and mi : m,2 = 11 : 9 is almost identical to that of a/m — 1 and mi : m.2 = 3:1 (cf. Figs. l9l and 1 10(1 . perhaps pointing 
to another relation in the evolution equations that could be used to reduce further the dimensionality of the total 
search space, analogous to the effective spins described in 0. 



VII. CONCLUSIONS 



We have derived a set of orbit-averaged post-Newtonian equations of motion for two spinning point masses through- 
out the binary inspiral. This makes it possible to model the evolution of the relative spin orientations over a large 
range of time scales. One important result is the discovery of nontrivial equilibrium solutions that allow the system 
to remain locked in a given orientation, exhibiting stable resonance behavior. Through the loss of energy and angular 
momentum to gravitational radiation, a binary system can undergo a transition from an unlocked to locked position, 
in turn significantly affecting the distribution function of the spins in parameter space. 

This behavior appears to be due in large part to the constraints placed on the orbital angular momentum and spin 
vectors within the precession framework. Conservation of total angular momentum J = L + S restrict the available 
regions of parameter space accessible to the system over short time scales. In this context, it may be more helpful 
to think of the spin interactions as "spin-orbit-spin" coupling. The precession of Si will change the orientation of 
L, which in turn affects the precession of S2. In an effort to better understand the origin of this behavior, we have 
repeated many of the evolution calculations, now with the direct spin-spin precession terms (order 1.5PN) turned off. 
This can be done trivially by removing the terms proportional to spin in equations (|2.10al2.10bll . Interestingly, the 
effect is almost negligible, giving qualitatively the same results shown in Figures [S] and El However, we should be 
careful to note that this is not quite the same condition as applied in 0, Q , where the analysis is limited to spin-orbit 
interactions simply by disregarding the spin of the smaller body. The approach here is more closely analogous to the 
effective spin method of Buonanno et al. • We believe their high level of success in using a single pseudo-physical 
spin may be closely related to the resonances described in this paper. 

Coupled with an estimate for the initial spin-orbit orientation from binary evolution codes, the analysis presented 
in this paper can be used to predict the distribution of spins relevant for gravitational wave detectors and waveform 
templates. By limiting the template library to the family of equilibrium solutions, we could greatly reduce the size of 
the parameter search space, in turn increasing the chances of signal detection. The spin orientation at small orbital 
separations, shortly before merger, could also be very useful information for studying the highly relativistic plunge 
regime, where the post-Newtonian approximation breaks down. Among other difficulties, the problem of black holes 
merging in numerical relativity is plagued by a lack of insight into the appropriate initial conditions. The results of 
this paper may serve to provide some of that insight. 

Directions for future work include modeling a vastly expanded parameter space and in turn studying the resonance 
behavior over a greater region of astrophysical interest. These studies should include both LIGO and LISA type 
sources, as well as continuing to investigate other evidence for spin-orbit and spin-spin interactions, such as the 
electromagnetic signatures of supermassive black hole mergers and binary pulsar systems. Similarly, comprehensive 
calculations of compact binary formation scenarios are needed to provide reasonable initial data for the inspiral 
evolution. 

As we have mentioned repeatedly throughout this paper, a major part of the detection strategy for gravitational wave 
observatories is that of template matching. We propose using the equilibrium solutions as a subset of the parameter 
space as a way to maximize the computational efficiency of matched filtering. A next step would be to quantify this 
presumption by calculating the fitting factors of waveforms from off-resonance and near-resonance systems compared 
to those of the exact equilibrium solutions, similar to the analysis of 0. This approach would also provide concrete 
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predictions for detection rates as a function of the source distributions and could be used to give stronger upper limits 
on binary coalescence rates. By using this set of physical templates, the intrinsic binary parameters of the compact 
system should be determined with greater precession than attainable by using a psuedo-physical template family. 

Aside from the potentially large impact on gravitational wave detection, the spin locking result is also of considerable 
interest for general dynamical systems exhibiting resonance behavior. We have made some progress in understanding 
the geometry of this resonance and how systems can be captured into it via radiation reaction. A more comprehensive 
study of the dynamics, including the possible identification of other resonances not yet discovered, could in turn 
predict new signatures to look for in gravitational wave signals. Of particular interest would be the possibility 
of finding overlapping resonance regions, classically an important breeding ground for chaotic systems |25|. One 
promising approach to the resonance problem is the use of conservative Hamiltonian equations of motion (e. g. 
[2ti l27l I2q 1. as opposed to the Lagrangian methods used here and in most other post-Newtonian calculations. The 
Hamiltonian approach more closely resembles the classical formulations of dynamic systems like the Solar System and 
coupled oscillators that show many examples of resonance behavior |24j , and thus could give us greater insight into 
the dynamics of compact binaries. 
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APPENDIX A: DERIVATION OF THE EVOLUTION EQUATIONS 



Assuming elliptical orbits with semimajor axis a and eccentricity e, we define a coordinate system with the origin 
at the occupied focus, e z = Ljv/|Ljv|, and aligned with the ellipse's pericenter. Then r is given by 



r = [r cos /, r sin /, 0] , 



(Al) 



where / is the true anomaly. From the standard binary relations for elliptical orbits, we define the specific angular 
momentum of the reduced two-body system as 



with 



I = r 2 f = \J ma(l — e 2 ), 
a(l-e 2 ) 



1 + e cos / 

and the total mass m = mi + mg. We can then change variables from dt — > df to get 

[a(l-e 2 )] 3 / 2 



dt 



,df. 



m}f 2 (l + e cos /) 2 

Averaging the r-dependent part of equations (|2.2al I2.2b(l over one binary period P gives 



1 



dt 



p I 3(*"S)* = 



1 



27r df(l + ecos/) . 

(b cos / + S u sin /)(cos fe x + sin fe v 



27to 3 / 2 J [o(l - e 2 )] 3 / 2 
S - (e z ■ S)e z 



2a 3 (l-e 2 ) 3 / 2 ' 
Then the orbit-averaged precession vectors are 



dtSlx = 



l 3(l_ e 2)3/2 



3 m 2 3 S2 • L 
2m7 ~ 2~L 



N 



Lat + -S 2 



(A2) 



(A3) 



(A4) 



(A5) 



(A6) 



and 
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As explained in the text, we are primarily interested in the relative orientation of the three vectors Ljv, Si, and 
S2, a problem that can be reduced to only four variables (L, 61,62, A0) in the case of circular orbits. To avoid certain 
sign ambiguities, it is actually more convenient to use five variables, defined as following (normalizing to m = 1): 

L 

z\ 

zi 
P 



1^1= /xr 1 / 2 , 


(A8a) 


Si • Jjpf 

C ° S6l= Sl L ' 


(A8b) 


S2 • 

cos 02 = , 


(A8c) 


a Si • S 2 

COS 6\2 = Q c , 


(A8d) 


L N • (Si x S 2 ) 

LS1S2 


(A8e) 



Under radiation reaction, the separation r will monotonically decrease throughout the evolution, allowing us to define 
a new time variable r = ro — r so that, from (|2.5I) . 



64/^ d 
5r 3 c?r 



In these coordinates, with an over-dot ' representing d/dr, the equations of motion are: 

,2 



L = 



/r 
2i' 
15 



z\ = 



Z-2 



= 



a 



( 1 giji 

128m" 2 



15 
128m 
15 

128m 
1 



L 

S2Z2 1 
L mi 



T m 2 mi S1Z1 - ^2^2 
aL I 1 

mi m,2 L 



ziz 2 (3 + f3(ziZ2 + z 2 zi) - z\z\ - z 2 z 2 - P(3 



From these functions, the original variable Acf> can be easily restored through 

Acj) = tan _1 [o;, (/3 — zi^)]. 



(A9) 

(AlOa) 
(AlOb) 

(AlOc) 
(AlOd) 
(AlOe) 

(All) 



For reference, equations IjAlOb -d) can be compared with similar results derived in 0, equations (2a-c) and 0, 
equations (17-19). 
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FIG. 1: Schematic diagram of the spin and orbital angular momentum vectors. The coordinate system is defined such that 
Lat is along the z-axis and (61, 62) are the respective angles between Ljv and (Si, S2). The projection of Si onto the x-y plane 
is defined to be along the x-axis so the azimuthal spin angles are 4>i — and 4>2 = A<f>. 
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FIG. 2: Equilibrium solutions of equation (13.511 for the case of maximally spinning black holes with normalized masses 
(mi = 0.55, rr%2 = 0.45) and spin alignment A(j> — 0°. The solution curves are labeled by their different values of the orbital 
angular momentum Ln, so that equilibrium systems evolving under radiation reaction move along subsequent curves with 
decreasing Ln (corresponding to r/m ~ 260, 150, 65, 16, 4). The labels 'A' and 'B' refer to the initial conditions for Figures^ 
andQJ), both with L N (t ) = 4. 
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FIG. 3: Equilibrium solutions of equation (13.511 for the case of maximally spinning black holes with normalized masses 
(mi = 0.55, rri2 = 0.45) and spin anti-alignment A<f) — 180°. The solution curves are labeled by their different values of the 
orbital angular momentum Ljv, so that equilibrium systems evolving under radiation reaction move along subsequent curves 
with decreasing Lm (corresponding to r/m ~ 260, 150,65, 16,4). Solid curves correspond to stable equilibrium orientations, 
while the dashed curves are the quasi-stable solutions described in the text. The label 'C refers to the unstable initial conditions 
of Figure with Ljv(to) = 1. 
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FIG. 4: Integrated phase difference [A^>(t)] to t between two precessing spin vectors. Examples are shown for initial conditions 
(a) near a stable equilibrium with 9i = 50°, 62 = 120°, Acj> = 0°, and Lm — 4; (b) far away from equilibrium with 9i = 
120°, 6 2 = 50°, A<f> = 0°, and L N = 4; (c) near a quasi-stable equilibrium with 9\ = 30°,6» 2 = 133°, A(j> = 180°, and L N = 1 
(compare with Figures [5] and |SJ . For all three cases, radiation reaction has been turned off. 
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FIG. 5: Sinusoidal projections of spin parameter space [612 = cos - (Si • S2),A^£>]. A sinusoidal projection maps spherical 
coordinates (8,<j>) to flat Cartesian coordinates ((f>sin8,9) with equal solid angles mapping to equal areas. The snapshots 
show a random sample of initial conditions with r(to) = 1000m evolving under radiation reaction at times when r/m — 
1000, 250, 125, 10. The masses of the two compact objects are similar, with mi = 0.55 and m2 = 0.45, and both have maximal 
spins. The initial spin of the larger mass is closely aligned with the orbital angular momentum: 0i(to) = 10°. 




FIG. 6: Sinusoidal projections of spin parameter space [612 = cos -1 (Si ■ S2), A0], as in Figure^] The snapshots show a random 
sample of initial conditions with r(to) = 1000m evolving under radiation reaction at times when r/m = 1000, 250, 125, 10. The 
masses of the two compact objects are similar, with mi = 0.55 and m,2 = 0.45, and both have maximal spins. The initial spin 
of the larger mass is misaligned with the orbital angular momentum: 61 (to) = 170°. 
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FIG. 7: Probability distribution for spin orientations (a) cos #12 and (b) A(j> near the end of inspiral [r/m = 10). The 
distributions are plotted for different values of 9\(to) with r(0) = 1000m in all cases. Both black holes are maximally spinning 
with similar masses mi = 0.55, m,2 = 0.45. When the larger spin is initially aligned with the orbital angular momentum 
(6*i < 90°), the final spins tend to align (A<j> — > 0°, cos #12 — ► 1). 




FIG. 8: Probability distribution for spin orientations (a) cos #12 and (b) A<f> near the end of inspiral. All parameters are as in 
Figure |7| yet with the larger spin initially misaligned with the orbital angular momentum [8± (to) > 90°]. In these cases, the 
final spins also tend to be anti-aligned with each other (A<f> — > ±180°, cos #12 — * —1), although cos #12 is less sensitive to the 
initial conditions, as can be seen from Figure |S| 
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FIG. 9: Probability distribution for spin orientations (a) cos #12 and (b) A(j> near the end of inspiral. The initial spin 
orientation has 8\ — 10° and uniform distribution of 82- Each system has maximal spins but a different mass ratio mi : 7712. 
The spin-locking resonance is clearly more important for nearly equal mass systems, as seen from equation 15,21 . 




FIG. 10: Probability distribution for spin orientations (a) cos #12 and (b) A(j> near the end of inspiral. The initial spin 
orientation has 81 — 10° and uniform distribution of 82. Each system has mi = 0.55, 7712 = 0.45 but different values for the 
black hole spin parameter a/m. The evolution is largely independent of the spin magnitude for moderately values of a/m > 0.5, 
while even systems with small spin parameters exhibit significant spin locking. 




